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Abstract. We study SU(2) lattice gauge theory with two flavors of Wilson fermion at non-zero chemical 
potential /i and low temperature on a 8^ x 16 system. We identify three regimes along the /i-axis. For 
/i < m^/2 the system remains in the vacuum phase, and all physical observables considered remain 
essentially unchanged. The intermediate regime is characterised by a non-zero diquark condensate and an 
associated increase in the baryon density, consistent with what is expected for Bose-Einstein condensation 
of tightly bound diquarks. We also observe screening of the static quark potential here. In the high-density 
deconfined regime we find a non-zero Polyakov loop and a strong modification of the gluon propagator, 
including significant screening in the magnetic sector in the static limit, which must have a non-perturbative 
origin. The behaviour of thermodynamic observables and the superfiuid order parameter are consistent with 
a Fermi surface disrupted by a BCS diquark condensate. The energy per baryon as a function of /i exhibits 
a minimum in the deconfined regime, implying that macroscopic objects such as stars formed in this theory 
are largely composed of quark matter. 



1 Introduction 

At large quark chemical potential /i, QCD is expected 
to undergo a transition from a confined nuclear matter 
phase to a deconfined quark matter phase, where the rel- 
evant degrees of freedom are quarks and gluons. It is now 
also generally believed that the quark matter phase at 
low temperature T is characterised by diquark condensa- 
tion: pairing of quarks near the Fermi surface gives rise to 
a number of color superconducting phases [T]P] |^ ]| ^ . The 
phase structure depends critically on the precise values 
of the diquark gap parameters and the effective strange 
quark mass at the relevant densities, and in the absence of 
a first-principles non-perturbative determination of these 
quantities, our knowledge of this region of the phase dia- 
gram will remain unsatisfactory. 

Lattice QCD is at present unable to address these 
problems directly, since the fermion determinant is no 
longer positive definite once n y^ 0, and cannot be used as 
a probability weight in the functional integral. There has 
been much progress in recent years in developing methods 
for the region of high T, low fj,, where the problem is less 
severe [H] . These methods are being extended in the direc- 
tion of higher n and lower T, but it is not clear at what 
point they break down, beyond which only unphysical re- 
sults will be obtained. 

Another approach is to study QCD-like theories where 
the fermion determinant remains real and positive even 
when /i 7^ 0. These can be used as a laboratory for investi- 
gating gauge theories at high density and low temperature. 
One such theory is 2-color QCD (QC2D) with an even 



number of flavors ^6 , where the quarks and antiquarks 
live in equivalent representations of the color group and 
can be related by an anti-unitary symmetry (the Pauli- 
Giirsey symmetry). This theory has been studied on the 

At /z = 0, the Pauli-Giirsey symmetry implies an ex- 
act symmetry between mesons and diquarks, which are 
the baryons of the theory. In particular, chiral multiplets 
will contain both mesons and baryons. For Nf = 2 for 
example, the pseudo-Goldstone multiplet consists of the 
pion isotriplct plus a scalar isoscalar diquark and antidi- 
quark. The diquark baryons can be expected to condense 
when fi > m.^/2, forming a superfluid ground state. In this 
respect, the theory is radically different from real QCD, 
where no gauge invariant diquark operator exists and the 
ground state at high density is superconducting. The na- 
ture of the superfluid ground state is however an inter- 
esting issue in its own right. For instance, an alternative 
approach to a superfluid order parameter in terms of an 
orthodox BCS description of diquark pairing at the Fermi 
surface has been given in jl4| . 

In the gluon sector, the differences between SU(2) and 
SU(3) are expected to be less important, and QC2D is a 
good setting for ab initio studies of gluodynamics in the 
presence of a background baryon density. Of particular 
interest is the issue of deconfinement at high density. Sig- 
nals of deconfinement were observed in simulations with 
Wilson PJ and staggered fermions TT , where correla- 
tions were found between the Polyakov loop and chiral 
or baryonic observables. However, the phase structure has 
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not been investigated in any further detail, and it remains 
unclear whether there is a confined nuclear matter phase 
with non-zero baryon number (as in QCD), or just a single 
phase transition. This will be investigated in the present 
paper. 



Most of the lattice studies so far have been conducted 
using staggered fermions, which have several potential the- 
oretical problems. Firstly, the pattern of global symmetry 
breaking is different to that of the continuum, so two-color 
staggered lattice QCD has a different Goldstone spectrum 
to continuum QC2D |S1|H1- Secondly, one staggered fla- 
vor corresponds to four continuum flavors, which is un- 
comfortably close to the Banks-Zaks threshold N?^ = 
136/49 « 5.6 where the second term in the /3-function 
changes sign, leading to a nontrivial infrared fixed point 
and absence of confinement and chiral symmetry breaking. 
Note that this prediction for N?^ is inherently perturba- 



tive, and may be unreliable since the fixed-point coupling 
is large for Nf ~ N?^ . In order to describe a single contin- 
uum flavor one usually takes the fourth root of the fcrmion 
determinant. It is not clear whether this procedure yields 
something corresponding to a local action of a single-flavor 
fermion field |15|. Even if it can be shown to be a valid 
procedure at /Lt = 0, obstacles remain for /x 7^ which may 
invalidate it [TC| . 



These problems are absent if Wilson fermions are used. 
We note that the Wilson formulation still admits a U(l)s 
global symmetry implying a conserved baryon number, 
and so a superfluid order parameter remains well-defined. 
In the chiral limit k — + Kc, the lattice Dirac operator's 
eigenvalue spectrum lies in the same chiral orthogonal en- 
semble universality class as the continuum theory. More- 
over, the consequences of explicitly broken chiral symme- 
try should be less severe, since they manifest themselves at 
the bottom of the Fermi sea, and hence become physically 
irrelevant at large /i. Because the Wilson formulation has 
the whole of the first Brillouin zone available to describe 
a single physical flavor, saturation artefacts set in at a 
larger value of /i than is the case for staggered ^^l- On 
the other hand one has to contend with a higher compu- 
tational cost, and for this reason, only a few studies using 
Wilson fermions have been performed to date |11|. 



In this article we will present results from a study of 
QC2D with two flavors of Wilson fermion at zero tem- 
perature and non-zero chemical potential. In section|21we 
present the lattice action and expressions for the principal 
bulk observables. The simulation parameters are given in 
section O along with results for the lattice spacing and 
pion and rho meson masses from simulations at /i = 0. 
The main results are given in section ^ Finally, in sec- 
tion we discuss the significance of our results and the 
prospects for further work. 



2 Lattice formulation and simulation 
2.1 Action and algorithm 

The Nf = 2 fermion action is given by 

S = 4JiM{^i)i>i + i^2M{fi)i;2 

-j7/ii(C75)T-2V5r + J7^r(C^75)r2V^i, (1 
where M(fi) is the usual Wilson fcrmion matrix 

V 



(2) 



The diquark source terms J, J serve a double purpose in 
lifting the low-lying eigenmodes in the superfluid phase, 
thus making the simulation numerically tractable, and en- 
abling us to study diquark condensation without any "par- 
tial quenching" . In principle the results should at the end 
be extrapolated to the "physical" limit J = J = 0. We 
will also introduce the rescaled source strength j = J/k. 
The fermion matrix has the following symmetries: 

75M(^)75 = Aft(-/,) (3) 

KM{n)K-^ = M*{fi) with K = Cj5T2, (4) 

where we have used the property T2U^T2 = U*. 

The last equation is the Pauli-Giirsey symmetry. This 
symmetry implies that det M(/i) is real, but not necessar- 
ily positive. However, with the change of variables ((> = 
—ip^CT2,(t) = C~^T2'ip2^ we can rewrite the action as 



«-(^^)("i^^/frt 



'^ 1 = <1'M<1'. 



(5) 



Hence positivity of detAl requires the product J J to be 
real and positive, which translates into the requirement 
that the diquark source term be antihcrmitian ^18j . 
Now use (El to write 



M^M 



Aft(^)M(/.)- 



Aft(-^)Af(-/i) + |J|2 



(6) 



The off-diagonal terms can be shown to vanish if J = J* 
using |j2Il; moreover the same identity applied to the lower 
block yields 



det M^M = [det{M\n)M{n) + JJ)f 



(7) 



It is therefore possible to take the square root analyti- 
cally, by using pseudofermion fields with weight (Af 'I'M -|- 
|Jp)~^. This has the advantage of (a) requiring matrix 
multiplications of half the dimensionality, and (b) permit- 
ting a hamiltonian evaluation and hence the use of an 
exact HMC algorithm. It is equivalent to the even/odd 
partitioning step used for staggered fermion gauge the- 
ories, but is more transparent since all lattice sites are 
physically equivalent, making the force term easier to im- 
plement. The trick was used in TF", though because the 
staggered version still requires a Pfafhan rather than a 
determinant, an HMD algorithm was used in that case. 
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2.2 Observables 

The quark number density is given by the timehke com- 
ponent of the conserved vector current: 

i 

+ M^)ho + l)e-''U}{x ~ i)i;,{x - t)) . (8) 



terms of the conformal anomaly 



The quark energy density can be defined in terms of a 
local bilinear very similar to ((SJ: 

e, = ^ K(i>,{x){"fo - l)e''Utix)->Jj,{x + i) 

■i 

- ^^{x)(-fo + l)e"^C//(x - i>»(a; - t)) . (9) 

Unlike Uq, however, this quantity requires both additive 
and multiplicative renormalisation as a result of quantum 
corrections. First, the vacuum contribution e^o must be 
subtracted. This correction can be obtained in the zero 
temperature thermodynamic limit from the relation 



1 



D 



e^ = - Trl-(^^)^=o 



(10) 



valid in _D-dimensional spacetime, or more directly, as we 
do here, by subtracting the measured value of £g(/x = 0) ~ 
0.3982(8). We have verified that this gives the same re- 
sult as using eq. (fTn| . The multiplicative correction re- 
sults from the renormalisation of the anisotropy factor 
^ — at /as under quantum corrections; it affects the quan- 
tity defined in ^ even in the isotropic limit ^ = 1, and 
must be determined either perturbatively or preferably 
non-perturbatively, requiring simulations with ^ 7^ 1. To 
our knowledge the perturbative correction has yet to be 
calculated for Wilson fermions; for staggered fermions it 
has been computed to be 1 — 0.1599C2(A^c)5^ HO], where 
C2 is the quadratic Casimir. This suggests that our results 
for £q should be rescaled downwards. Since the correction 
is independent of /i, however, this is only an overall nor- 
malisation factor. 

With the standard Wilson lattice gauge action em- 
ployed in this study the gluon energy density may be de- 
fined as a local observable 



1 



mNt 



at 



dSq 

dat 



^^TT{nt-Os), (11) 



where Dj, D^ are timelike and spacelike plaquettes respec- 
tively. Once again, this requires renormalisation due to 
quantum corrections; the dominant correction factor has 
been calculated in perturbation theory to be 1— 0.1466C2g^ 
|?T|. suggesting that bare Eg data should again be rescaled 
downwards by a /i-independent factor. 

The final thermodynamic observable we can discuss is 
the trace of the energy-momentum tensor, expressible in 



3p 
1 



Nmt 



da 



da 



dlnZ 

LCP^P 

d\nZ 



dn 
da 



dh-iZ 



LCP 



Ok 



(12) 



LCP 



dj 



where it is understood that beta-functions are evaluated at 
/I = T = along lines of constant physics, so dimension- 
less ratios of physical quantities are cutoff-independent 
(the derivation of this equation for /i ^^ is sketched in 
|22|). Baryon number symmetry implies that limj_>o dj /da ■ 
and hence the third term can be neglected. The local ob- 
servables required at /i 7^ are then 

1 dlnZ 3f3 ^ , 

The final bulk observables we compute are the diquark 
condensate. 



(9g) = -(^iif^r-7^r^vi), 



(15) 



which is an order parameter for the vacuum-to-superfluid 
transition; and the Polyakov loop N~^{TrL), which in 
pure gauge theories is an order parameter for the decon- 
finement transition. 



3 Simulation parameters and physical scales 

Exploratory simulations at zero chemical potential were 
performed at a range of different values for (3 and k ^Sl > 
but for this study the parameters P = 1.7, k = 0.178 were 
selected. All the simulations have been performed on a 
8^^ X 16 lattice. The full set of parameters used is given in 
tabled Configurations were saved every 4 trajectories. 

In order to determine the scale, we computed pi and 
rho meson correlators and the static quark potential. The 
mesons were computed using point sources (no smearing), 
and as fig. ^ and table [5] clearly show, on our relatively 
small lattice this did not allow a precise determination 
of the mass. In a future in-depth study of the hadron 
spectrum it will be desirable to use smeared sources and 
variational techniques, but at this point we are primarily 
interested in a rough idea of the hadronic scales. These 
numbers indicate that we may expect an onset transition 
at fio ~ '7i7r/2 ~ 0.4a~^, but since there is no separa- 
tion of scales between the Goldstone (pion) and the non- 
Goldstone (rho), we should not expect chiral perturbation 
theory to be quantitatively valid at any point. 

The static quark potential was computed using APE 
smeared Wilson loops with spatial separations near the 
diagonal to minimise lattice artefacts. The minimum time 
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Ntr 



{£) 



dt 



iVc, 



Range 



X'/Ndf 



xVN^f 









2000 


1.0 


0.0125 


385 


85% 


0.1 





608 


1.0 


0.0125 




82% 


0.2 





560 


1.0 


0.01 


785 


75% 


0.25 


0.04 


266 


0.5 


0.0075 


740 


86% 


0.3 





600 


1.0 


0.004 


1350 


80% 




0.02 


188 


0.5 


0.006 


1170 


80% 




0.04 


190 


0.5 


0.0075 


965 


75% 




0.06 


276 


0.5 


0.0075 


775 


86% 


0.35 





500 


0.1 


0.0005 


1650 


90% 




0.02 


400 


1.0 


0.004 


1615 


75% 




0.04 


500 


1.0 


0.005 


1090 


85% 


0.4 


0.04 


148 


1.0 


0.005 


1165 


80% 






302 


0.5 


0.006 


1235 


77% 


0.45 


0.04 


252 


0.5 


0.0042 


1275 


82% 






64 


1.0 


0.004 


1200 


88% 


0.5 


0.02 


554 


0.5 


0.003 


2565 


78% 




0.04 


44 


1.0 


0.004 


1240 


74% 






300 


0.5 


0.0045 


1270 


77% 




0.06 


304 


0.5 


0.005 


900 


86% 


0.55 


0.04 


308 


0.5 


0.0038 


1340 


83% 


0.6 


0.04 


112 


1.0 


0.0033 


1290 


83% 






283 


0.5 


0.004 


1375 


80% 


0.65 


0.04 


276 


1.0 


0.003 


1310 


85% 


0.7 


0.02 


400 


1.0 


0.002 


2525 


80% 






440 


0.5 


0.0025 


2760 


75% 




0.04 


136 


1.0 


0.003 


1330 


85% 






530 


0.5 


0.0035 


1450 


78% 




0.06 


280 


0.5 


0.0037 


970 


85% 


0.75 


0.04 


180 


1.0 


0.003 


1390 


80% 


0.8 


0.04 


292 


1.0 


0.0028 


1410 


80% 






292 


0.5 


0.003 


1520 


74% 


0.9 


0.04 


264 


0.5 


0.003 


1610 


72% 


1.0 


0.04 


272 


0.5 


0.0028 


1660 


75% 


1.3 


0.04 


249 


0.5 


0.0033 


1070 


88% 


1.4 


0.04 


404 


0.5 


0.005 


620 


82% 


1.5 


0.04 


400 


0.5 


0.01 


290 


80% 


1.6 


0.04 


400 


0.5 


0.01 


215 


90% 


1.75 


0.04 


250 


0.5 


0.01 


200 


94% 


2.5 


0.04 


260 


0.5 


0.012 


260 


97% 


Tabk 


! 1. Simulation 


parameters. {£) denotes the 


average 


jectory length, dt is 


the molecular dynamics timestep, 



is the number of conjugate gradient iterations needed in the 
molecular dynamics evolution, and ace is the acceptance rate. 




3-4 


0.809(2) 


0.58 


0.882(4) 


0.71 


0.917(3) 


3-5 


0.800(2) 


0.97 


0.870(3) 


1.4 


0.919(3) 


3-6 


0.795(2) 


1.3 


0.860(3) 


2.5 


0.924(3) 


4-5 


0.789(2) 


0.81 


0.854(4) 


1.7 


0.924(3) 


4-6 


0.786(2) 


0.63 


0.846(3) 


1.9 


0.929(3) 


5-6 


0.783(2) 


0.80 


0.836(4) 


3.0 


0.937(3) 



Table 2. Fit ranges and fitted masses for pi and rho mesons. 



> 




R/a 

Fig. 2. Static quark potential for fj, = j — 0, together with a 
fit to the Cornell potential. 



separation used was 2. Figure |21 shows the potential to- 
gether with a fit to the Cornell potential V{R) = aR + 
e/R + C up to oiJmax — 4.0. We find for the string tension 
cra^ = 0.218(8), which gives a 
spacing. 



0.223(4) fm for the lattice 



4 Results at yu 7^ 

4.1 Model Considerations 

Our expectation as fi is increased from zero at T ~ is 
that the system will remain in the vacuum phase until an 
onset occurs at /^ = /io — Mught/A^c, where Mught is the 
mass of the lightest baryon in the physical spectrum. In 
QC2D this lightest baryon is a scalar diquark state in the 
same chiral multiplet as the pion: hence fig = W7r/2, and 
in the chiral limit the onset transition is well described 
by an effective chiral model (xPT) in which only pseudo- 
Goldstone pion and diquark degrees of freedom are re- 
tained. A mean-field treatment of xPT has yielded quan- 
titative predictions for chiral and superfluid condensates, 
quark number density, and Goldstone spectrum as /x is 
increased beyond onset 23 . Our starting point for ther- 
modynamics is the result for quark number density: 



rixPT 



fo, 

|8A^/FV( 



1 



M > Mo- 



(16) 



Fig. 1. Pi and rho effective masses for n — j = 0. 



Here F is the pion decay constant, a parameter of the 
xPT which can be extracted in principle from the pion 
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correlator measured at /i = 0. The pressure follows im- 
mediately from integration of the fundamental relation 
Uq = dp/dfi\T,V- 






(17) 



Now, if f2{fJ-, T, V) is the thermodynamic grand potential, 
then in the grand canonical ensemble f2 = —pV. In the 
T = limit we can therefore extract the energy density 
via f2/V = Eq — ^Uq, implying 



e^PT = 4A^/F^ ^ 



9 



Sa^o 



(18) 



We can also extract the trace of the energy-momentum 
tensor 5 = Sq — 3p: 






(19) 



and note that 5^pt is positive for Ho < fJ- < \/ofJ-o- 

The result IjlSI) for e^pT coincides with that derived 
in p3" using the assumption that diquarks Bose-condense 
in the ground state, remain degenerate with the pion for 
jj, > ^o, and that their self- interactions are weak due to 
their Goldstone nature. One infers Sq ~ ^mT^nq = ^ionq. 
A more refined treatment taking account of corrections to 
the dilute ideal Bose gas from the weak repulsive interac- 
tion yields PH] 



Mo"-g 



64F2iV 



/ 



(20) 



where the dots denote corrections from subleading terms 
in the chiral expansion. Using /i = dsq/dnq, this can be 
seen to be consistent with (|16f) linearised about fi = iio'- 



l^PT 



32NfF\ti~tioy, 



(21) 



the same approximation predicts p^pt to vanish as (/i — 
/Zo)^, consistent with a second order transition in the Ehren- 
fest scheme. 

Next we turn to the deconfined phase expected at large 
fi. The obvious starting point is the Stefan-Boltzmann 
prediction for the number density of massless quarks: 



nsB 



NfN, 
37r2 



-M 



(22) 



This is obtained simply by populating a Fermi sphere 
of radius kp = /z, with every momentum state occupied 
by 2NfNc non-interacting quarks. Other thermodynamic 
quantities follow immediately: 



£SB = 3p, 



SB 






SSB = 0. 



(23) 



An estimate for the chemical potential fid at which de- 
confinement takes place can be obtained by equating the 
free grand potential densities, or equivalently the pressures 
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2.5 
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Fig. 3. Ratio of thermodynamic observables from xPT and 
free quarks 



given in (|17|l and (|23|l . Since p^pt > Psb for /i < /i^ is 
required for thermodynamic stability, we find /i^ given by 
the largest positive real root of 



^^d 



' ^^F\hriA-^^l) = ^■ 



(24) 



This estimate takes no account of any non-Goldstone states 
in the hadron spectrum, or of any gluon degrees of free- 
dom released at deconfinement. In fig. |21we plot the ra- 
tios n^px/nsB, ^xPt/^sb and PxPt/psb as functions 
of n/no for the choice F^ = Nc/Gir'^, corresponding to 
fid ~ 2.288/io. Since n^pp < nsp at this point, this 
naively simple approach predicts a first order deconfining 
transition — note also that S^pp < at deconfinement. 
Also shown is the ratio of the speed of sound 



'JxPT 



1 



(25) 



to the Stefan-Boltzmann value vsb = l/v3. 
4.2 Thermodynamics results 



The raw data for bosonic observables (spatial and tempo- 
ral plaquettes, Polyakov line) are tabulated in tabled and 
for fermionic bihnears {{'iij'ip) , Uq ^, Sq © and (qq) (ESJ) 
in table^l All quark observables are normalised to N, col- 
ors and Nf flavors, with Nf = N, = 2. In this section we 
outline the analysis needed to extract bulk thermodynam- 
ical quantities and condensates. 

The most straightforward observable to analyse is the 
quark density n^, which as the timelike component of 
a conserved current requires no renormalisation due to 
quantum corrections. There may, however, still be lattice 
artefacts in both UV and IR regimes, as illustrated by the 
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j 


M 


Tr(a, + 


□t)/2iVc Tr(nt - D,)/iV, 


TtL/Nc 


0.00 


0.00 


0.4738(1] 


-0.0002(1) 


-0.0014(05) 




0.30 


0.4739(3] 


0.0001(2) 


-0.0008(14) 


0.02 


0.30 


0.4742(3] 


0.0002(2) 


-0.0012(15) 




0.35 


0.4755(4] 


-0.0005(3) 


0.0005(24) 




0.50 


0.4804(5] 


0.0024(2) 


0.0020(13) 




0.70 


0.4773(3] 


0.0089(2) 


0.0127(11) 


0.04 


0.25 


0.4735(3] 


0.0002(3) 


-0.0026(14) 




0.30 


0.4753(4] 


0.0007(4) 


-0.0013(20) 




0.35 


0.4764(4] 


0.0004(2) 


0.0003(13) 




0.40 


0.4784(4] 


0.0012(3) 


0.0020(14) 




0.45 


0.4792(3] 


0.0020(4) 


-0.0016(18) 




0.50 


0.4799(4] 


0.0025(3) 


0.0044(15) 




0.55 


0.4793(4] 


0.0037(3) 


0.0004(17) 




0.60 


0.4794(3] 


0.0051(3) 


0.0029(21) 




0.65 


0.4778(3] 


0.0069(2) 


0.0065(11) 




0.70 


0.4773(4] 


0.0092(3) 


0.0094(16) 




0.75 


0.4752(2] 


0.0122(2) 


0.0188(19) 




0.80 


0.4719(2] 


0.0172(3) 


0.0292(12) 




0.90 


0.4648(3] 


0.0299(3) 


0.0719(17) 




1.00 


0.4549(3] 


0.0442(3) 


0.1393(19) 




1.30 


0.4277(4] 


0.0479(3) 


0.4286(20) 




1.40 


0.4142(3] 


0.0259(4) 


0.3216(11) 




1.50 


0.4060(2] 


0.0128(3) 


0.1019(13) 




1.60 


0.4052(2] 


0.0114(3) 


0.0224(18) 




1.75 


0.4052(2] 


0.0105(4) 


0.0018(18) 




2.50 


0.4053(3] 


0.0115(4) 


0.0008(13) 


0.06 


0.30 


0.4767(3] 


0.0003(2) 


0.0004(16) 




0.50 


0.4791(4] 


0.0025(2) 


-0.0035(14) 




0.70 


0.4772(3] 


0.0095(4) 


0.0094(16) 



Table 3. Raw data for gluonic observables 



equivalent quantity for free fields on a N^ x Nt lattice: 



where 



ANfNc 



T.TT 



i sin fco [Li cos ki - ^] 



fc [^-E^cosfc^] -KX;^sin2fc^ 

(26) 



fco - i^i = It ("0 ' ^ ' 



27rnt^ 



Nt 



Ifi, 



i/= 1,2,3. 



(27) 
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For free massless quarks k = |- In the large-/x limit n^ssif^) 
saturates at a value 2NfNc per lattice site. This is an arte- 
fact of non-zero lattice spacing, which we will discuss in 
more detail in section 14.31 below. For T — the corre- 
sponding continuum relation is (f^ . Fig.0]plots n^sB/f^^ 
for several system sizes, and shows there are significant 
departures from the continuum result both at small fi as 
a result of finite Ng, and at /xa > 0(1) as a result of 
non-zero lattice spacing, species doubling etc. Of partic- 
ular interest are the pronounced wiggles seen especially 
on the 8^ x 64 curve. These arise due to departures from 
sphericity of the Fermi surface in a finite spatial volume, 
and are visible whenever the temperature is much smaller 
than the mode spacing, i.e. Nt 3> Ng |24| . 

The main implication in the window /i € (0.25,1.0) 
where our analysis is focussed is that n}gg lies systemati- 
cally above the continuum value, with the ratio increasing have a similar shape, with a maximum departure over the 



Fig. 4. n^gg/fi'^ for free massless Wilson fermions. 



for /i > 0.5. To correct for this lattice artefact, we have 
chosen to plot the ratio nq/ri^g^, shown as a function of 
fi in fig. The free field density has been calculated for 
massless quarks, i.e. with k = i. In order to assess the 
systematic error due to our lack of accurate knowledge of 
the non-zero quark mass, we repeated the free field calcu- 
lation with K — 0.120, corresponding to a conservatively 
large niga — 0.167, and found the resulting ".^^(/i) to 
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j 


M 


rig 


Eg 


(ilj-iP) 


(W> 


0.00 


0.00 


-0.0004(04) 


0.3980(4) 


14.409(2) 






0.30 


-0.0008(22) 


0.3994(22) 


14.399(7) 


- 


0.02 


0.30 


0.0012(12) 


0.4028(16) 


14.389(5) 


0.0219(2) 




0.35 


0.0030(24) 


0.4072(24) 


14.364(3) 


0.0293(3) 




0.50 


0.0382(24) 


0.4496(40) 


14.254(10) 


0.0571(4) 




0.70 


0.1506(32) 


0.5348(32) 


14.178(6) 


0.1099(5) 


0.04 


0.25 


0.0036(14) 


0.4012(16) 


14.396(8) 


0.0343(1) 




0.30 


0.0102(14) 


0.4126(16) 


14.360(4) 


0.0400(2) 




0.35 


0.0144(12) 


0.4188(18) 


14.332(4) 


0.0469(3) 




0.40 


0.0252(16) 


0.4358(16) 


14.292(8) 


0.0549(3) 




0.45 


0.0332(24) 


0.4422(30) 


14.272(8) 


0.0634(5) 




0.50 


0.0476(22) 


0.4558(22) 


14.244(8) 


0.0713(5) 




0.55 


0.0648(30) 


0.4702(30) 


14.236(8) 


0.0817(5) 




0.60 


0.0898(30) 


0.4902(30) 


14.196(8) 


0.0932(4) 




0.65 


0.1238(26) 


0.5146(28) 


14.180(4) 


0.1084(4) 




0.70 


0.1660(24) 


0.5476(26) 


14.156(4) 


0.1248(4) 




0.75 


0.2302(34) 


0.5972(32) 


14.136(8) 


0.1461(5) 




0.80 


0.3354(38) 


0.6850(32) 


14.084(8) 


0.1701(5) 




0.90 


0.6664(56) 


0.9716(60) 


13.916(8) 


0.2222(9) 




1.00 


1.2602(74) 


1.5168(72) 


13.560(8) 


0.2654(12) 




1.30 


5.220(6) 


5.284(6) 


12.036(8) 


0.0956(3) 




1.40 


7.040(4) 


7.060(4) 


11.312(4) 


0.0594(1) 




1.50 


7.914(4) 


7.916(4) 


10.944(4) 


0.0346(3) 




1.60 


7.996(3) 


7.996(3) 


10.916(4) 


0.0255(1) 




1.75 


8.000(4) 


8.000(4) 


10.916(4) 


0.0207(1) 




2.50 


8.000(2) 


8.000(2) 


10.924(4) 


0.0161(1) 


0.06 


0.30 


0.0130(14) 


0.4184(14) 


14.331(6) 


0.0554(2) 




0.50 


0.0530(22) 


0.4604(22) 


14.230(8) 


0.0855(8) 




0.70 


0.1780(30) 


0.5574(30) 


14.131(8) 


0.1397(10) 



Table 4. Raw data for quark observables (no zero subtraction applied to e. 
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Fig. 5. Ratio of thermodynamic observables to free field values 
versus /i for j — 0.04. Open symbols show extrapolations to 



fi range of interest of 0(20%). Where data are available 
we have plotted as open symbols the result of linearly ex- 
trapolating J ^ 0; at /i = 0.3 this yields Ug « 0, consistent 
with being below the onset expected at fio — "t-tt/S. For 
H = 0.5, 0.7 the extrapolation results in a downwards cor- 
rection of 0(20%). We also note that for /i > 0.5 the ratio 



is roughly constant and greater than one, which is plau- 
sible if Stefan-Boltzmann scaling sets in at large /i, but 
with a Fermi momentum kp > fi. Physically, this would 
result from degenerate quark matter with a positive bind- 
ing energy arising from interactions. Note that a non-zero 
quark mass has the opposite effect, tending to raise fi over 
kp. 

In fig.El we also plot the unrenormalised quark energy 
density £q/£^sB versus ^, where e'^^ is evaluated using a 
formula similar to H26|l and £5^* is given in (|23|l . Note 
that systematic errors due to non-zero quark mass and in- 
correct vacuum subtraction are potentially larger in this 
case, but have maximum impact at the lower end of the fj, 
range of interest. Once again, a, j ^ extrapolation has 
been done where possible, and the signal for /x = 0.3 is 
consistent with the vacuum. As for quark number density, 
the ratio £q/£^sB tends to a constant at large /x, which 
we interpret as evidence for the formation of a Fermi sur- 
face. This limit is approached from above, however, and 
in the range 0.35 < M ^ 0.5 Sq/esB actually peaks. The 
implications of this are discussed in sectional 

The pressure p is best calculated using the integral 
method [2H1: i-e. p{lj) = /^ nq{fj,')dfj,' . Note that since 
the only /^-dependence comes from the quark action, this 
expression gives in principle the pressure of the entire sys- 
tem including both quarks and giuons. In practice for data 
taken away from both continuum and thermodynamic lim- 
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Fig. 6. Pressure verus p. using two diflterent integral methods ^'S- '^ ' Comparison of bare quark and gluon energy densities 

versus p. 



its we should make some attempt to correct for artefacts: 
we have experimented with two different ad hoc proce- 
dures: 



V 



Zn,{f,')df/ 



J u 



71-; 



II 



-(m) 

PSB Jf,^ '-SB 

PSB 



■{fi')nq{n')dn' 



JfJ-a 



SB 



{^i')d^i' 



(28) 
(29) 



SB are defined in eqs. (|22I26|I respectively. 



where n™^*, 

Ultimately, data from different physical lattice spacings 
will be required to determine which method is preferred. 
Using an extended trapezoidal rule to evaluate the integral 
on the j = 0.04 data of fig. |S1 we have estimated both 
p^ {fj.) and p^^ {n) and plot the results in figEJ In both cases 
the pressure rises monotonically from near zero at onset, 
but for method II there is some suggestion of a plateau in 
P/PSB for /i > 0.6. By jjl ~ 1.0 the two methods agree, with 
p w 2psB- This is consistent with the ratio Uq/nsB in the 
same regime, again suggesting a Fermi surface with kp > 
/x. For comparison with the other fermionic observables 
we have plotted p^^ /psb in fig- HI 

Next we turn to the gluon sector. Since at the cur- 
rent lattice spacing the perturbative correction to the bare 
gluon energy density is likely to be inadequate, in fig. |7| 
we content ourselves with plotting unrenormalised data for 
both Eg/fj,'^ and for comparison Sq/ fi^, at j = 0.04. There 
is no evidence in our data for any significant variation 
with i for gluonic observables. The plot illustrates: (i) the 
significant impact of rescaling the fermionic data by e'^^j 
(cf. fig-EJ; (a) for all ^ > fio, the gluonic contribution is a 
significant fraction of the total energy density; (in) most 
strikingly. Eg scales as /i^ over the whole range of /i. On 
dimensional grounds, this is the only physically sensible 
possibility if thermal effects are negligible. We reiterate 
that gluonic contributions to thernrodynamic observables 
are not present in the SB formulae (|23|l . so that observa- 
tions (ii) and (Hi) are non-trivial predictions of the simu- 
lation. Note that our lack of knowledge about the energy 



density renormalisation factors prevents us from making 
any quantitative estimate of the gluon energy density as 
proportion of the total energy density. 

The quark and gluon contributions to the conformal 
anomaly p4ll3|l are plotted in fig. |H1 As for energy den- 
sity, the raw data requires a zero-/i subtraction. To extract 
S, we also need information on the lattice beta-functions, 
which require an extensive simulation canrpaign and are 
not yet available. For this reason we restrict ourselves to 
qualitative remarks. The quark contribution shown in the 
upper panel of fig. |S1 increases monotonically from /i = 0. 
The sign of dn/da is found to be negative in QCD simu- 
lations J26II27J (although it can change sign as the chiral 
limit K —f Kc is approached away from the continuum limit 
(26"). This is in accord with what we have found in our 
SU(2) simulations at various P,k ^]. We conclude that 
the quark contribution to d is negative. The gluon con- 
tribution, by contrast, starts positive for fJ-o '^ fJ- '^ 0.6 
before changing sign to decrease montonically for fi > 0.6 
(recall 8(3/ da < due to asymptotic freedom). 

In fact, we expect the gluon contribution to dominate, 
since 6k/ Sp « 0(0.01) along lines of constant physics j27j . 
This is clearly required for consistency with fig. |S1 since 
the sign of S must coincide with the sign of e/esB —p/psB , 
which looks positive for 0.4 < /i < 0.6. We also get a clue 
about the unknown renormalisation factor for e^, since if 
5 is to go negative at large /j,, eron/ess < p/psb in this 
regime. 

The non-monotonic behaviour of the plaquette has been 
predicted in a xPT study in which asymptotic freedom of 
the gauge coupling is taken into account |2H1, and has 
also been observed in recent simulations with staggered 
fermions ^21 ■ It can be understood in terms of Pauli block- 
ing. In quark matter as /i increases, the number oi qq pairs 
available for vacuum polarisation corrections to the gluon 
decreases, since only states close to the Fermi surface can 
participate. In the limit of complete saturation (i.e. one 
quark of each color, flavor and spin per lattice site) the 
gluon dynamics hence resembles that of the quenched the- 
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Fig. 8. Top: k~'^ Tr{l- (ipip)), andhottom: 3l3/NcTT{Dt + a^), 
as functions of /i at j = 0.04. 



ory, so Tr n(/i ^ oo) < Tr n(^ = 0). Assuming a smooth 
passage to the limit we deduce J < at large /.t. 



To summarise the thermodynamic information, we have 
been able to extract Uq (/z) and p(/Lt) directly from the sim- 
ulation — any remaining UV and IR artefacts in fig. [S] 
can be controlled by simulations closer to thermodynamic 
and continuum limits. The status of £q,£g is less secure, 
because these quantities require renormalisation by an as 
yet undetermined factor; similarly the trace of the energy- 
momentum tensor 6 requires knowledge of the lattice beta- 
functions. In each case, though, the rescaling factor is /i- 
independent, so the shape of the data in figs.[Sland|Slis cor- 
rect. However, the behaviour of e, in fig.|U and Sg in fig.|Hl 
give a strong hint of two qualitatively distinct high density 
regions: (i) for A^o ^ M ^ 0.65 Eg « 0{5esB), and Sg > 0; 
(ii) for fj, > 0.65 Eg/esB ~ nq/nsB, Sg < 0. All thermo- 
dynamic quantities seem to have the same /i-scaling as 
their Stcfan-Boltzmann counterparts in this higher den- 
sity regime, although Ug/nss ~ p/psb ~ 2. 



6- 



♦ ♦■♦■■ 



-- 8 xl6SB 
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Fig. 9. Quark number density n, and Polyakov loop L vs. n 
for values up to p = 2.5. All lattice data are at j = 0.04. 



4.3 The approach to saturation 



In figureElwe plot riq{fi) including values of /i up to 2.5, as 
tabulated in tables Q] 13 and ^ together with both contin- 
uum (|22|l and lattice 1)26(1 expectations for non-interacting 
quarks. On any system with a non-zero lattice spacing, 
there is a point at which every lattice site is occupied by 
the maximum value 2NcNf allowed by the Pauli exclu- 
sion principle. For free Wilson fermions this occurs for 
lia « 2.0; for the interacting quarks studied here the sat- 
uration threshold drops to /la « 1.5. We also see that the 



relation 



2"5B 



continues to hold all the way up to 



the threshold; there is no sign of asymptotic freedom. In 
this respect, the situation bears some similarity to that 
of QCD at temperatures between Tc and STc, where lat- 
tice simulations have also uncovered a deconfined, but still 
strongly interacting system. The nature of this system is 
quite different, however, as in the high-temperature case 
a slow approach towards Stefan-Boltzmann predictions is 
observed, while in the present case the strong binding en- 
ergy remains unchanged in the entire domain studied. 

As discussed above, in a saturated system, virtual qq 
pairs are suppressed, leading to the expectation that the 
gluodynamics should be that of the quenched theory. In 
fact, inspection of table O shows this is a little simplis- 
tic, since Eg remains non- vanishing even once saturation is 
complete. Strong coupling considerations suggest that the 
saturated system has a quark-induced effective action that 
can be expanded in even powers of the Polyakov loops. 
The resulting Scs distinguishes between □<, and □( hence 
yielding Eg ^ 0, but respects the global Z2 centre symme- 
try of the quenched action, consistent with (TrL) — 0. 

Just below saturation, Tr L rises rapidly. In this regime 
the theory resembles a p-type semiconductor, in that the 
low energy excitations are de-confined holes. The holes 
appear weakly interacting: (qq), which also measures the 
density of hole-hole pairs, is small, and there is little evi- 
dence from Ncg (table nj for a light bound state. 
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Fig. 10. The Polyakov loop L and the superfluid diquark con- 
densate {qq) / y? for j = 0.04 as function of /i. Open symbols 
show {qq) / tJ? extrapolated to j = 0. 
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Fig. 11. Superfluid condensate {qq) vs. j for various /i 



4.4 Order parameters 



Might it be possible to reconcile the behaviour reported 
in section 14.21 with the models of section 14. If To eluci- 
date this question we next review order parameters for 
superfluidity and deconfinement. In fig. 1101 we plot both 
the superfluid condensate (qq) given in (|15|l . rescaled by a 
factor /i^^, and the Polyakov loop L. Both show a marked 
change of behaviour at /i ~ 0.6; for fj, greater than this 
value {qq)/fJ''^ is approximately constant, while L rises 
from zero. To interpret the superfluid condensate, we first 
need to compare data taken at varying j, shown in fig. 1111 
together with a linear extrapolation to j = 0. The slight 
non-linearity of the n = 0.3 data suggests that below on- 
set, \imj^o{qq) = 0, whereas for fi — 0.5, 0.7 the data 
extrapolate to a non- vanishing intercept, implying a non- 
vanishing condensate and hence spontaneous breaking of 
baryon number symmetry, i.e. superfiuidity. Extrapolated 
values are plotted as open symbols in fig. ^| 




Now, in the Bose condensate phase expected for /Zo < 
H < fid, xPT predicts [2S| 



(30) 



which is a monotonically decreasing function for ^ > \/2^o- 
For a Fermi surface perturbed by a weak qq attractive 
force, by contrast, the order parameter counts the num- 
ber of Cooper pairs condensed in the ground state, which 
all originate in a shell of thickness A around the Fermi 
surface, where A is the superfluid gap: hence 



('Z'Z) 



BCS 



cx A. 



fi- 



(31) 



The data of fig. 1 101 support this scenario for /i > 0.6, with 
A independent of /z. Further support for the importance 
of quark degrees of freedom at large fi comes from the 
Polyakov loop, which rises smoothly from zero at /i ~ 0.65. 
We thus tentatively assign fi^ ~ 0.65, marking a transition 
from confined scalar "nuclear" matter to deconfined quark 
matter. In condensed matter physics parlance this tran- 
sition would be characterised as one from BEC to BCS. 
Exposing the detailed nature of the transition will require 
many more simulations using a variety of source strengths, 
lattice spacings, and spatial volumes. 



4.5 Static quark potential 

The screening effect of the dense medium can be further 
investigated by studying how the static quark potential 
changes as the chemical potential increases. In figure IT^ 
we show the static quark potential for various values of 
fi and j — 0.04. Up to /i = 0.3 there is no change from 
the fj, = potential, while in the intermediate phase we 
see clear evidence of screening due to a nonzero baryon 
density. This is even clearer in fig. 1131 where the /i = j = 
potential has been factored out. 

In the deconfined phase a new pattern emerges, where 
the short distance potential is strongly modified, while at 
long distances we appear to see an increase with increasing 
fi rather than a decrease as expected. There are however 
indications that the long-distance screening at fi — 0.7 
may become slightly stronger as j -^ 0, as fig. 1141 shows. 
No dependence on j is seen at /z = 0.3 or /i = 0.5. We do 
not have an understanding of these effects at present, al- 
though it is possible that lattice artefacts may contribute 
to the short-distance modifications. This can be resolved 
by going to finer lattices. Likewise, larger lattices will 
be required to determine the long distance potential to 
greater accuracy. 



4.6 Gluon propagator 

We have computed the gluon propagator by fixing the con- 
figurations to Landau gauge using an overrelaxation algo- 
rithm to a precision \d^A^\'^ < 10" ^°. The approach and 
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Fig. 12. Tlie static quark potential for various values of the 
chemical potential ^ and diquark source j — 0.04, together 
with the fi — j — potential shown in fig. |2| 




Fig. 15. Gluon propagator at /i = j — 0. 
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Fig. 13. The static quark potential for various values of the 
chemical potential ^ and diquark source j — 0.04, divided the 
fi = j — potential shown in fig. |5] 
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notation is analogous to that of 29 . At nonzero chemical 
potential, which defines a preferred rest frame, the gluon 
propagator can be decomposed as 



D^.Aq) = PntiQ)DM{q'o,q' 



+ P;^M)DE[qlq^) + PUq)DL{ql q' 



where 



qiqj 



M , 



jM 



^'(9) = (-5.,-^); P^iq)^P^=0; 



PuJq) 



PL{q) 



qu 



qtj,qv 
ql + (i 



q 

ql + (i 






PuUq); 



2 



1, 



2 



(32) 

(33) 
(34) 
(35) 
(36) 



Dm is the magnetic (spatially transverse) gluon propaga- 
tor and De is the electric (spatially longitudinal) prop- 
agator, while the longitudinal propagator D^ is zero in 
Landau gauge. 

Figure [Tsl shows the gluon propagator at /i = j = as 
a function of 4-momentum q^ = QQ+q"^. Since the Lorentz 
symmetry remains unbroken here there is only one form 
factor, and the small splitting between electric and mag- 
netic gluon is most likely a finite volume effect, similar to 
the asymmetric finite volume effects on the tensor struc- 
ture observed in the SU(3) gluon propagator [221 ■ 

Figure El shows the gluon propagator at /.t = 0.7, j = 
0.04 as function of the 4-momentum q^, compared to the 
vacuum gluon propagator. We find that the electric prop- 
agator remains virtually unchanged at this point, while 
some modifications can be seen in the magnetic prop- 
agator. In particular, the Lorentz or 0(4) symmetry is 
clearly broken since different values for q^ give different 
"branches" . 

In order to see the /i-cvolution more clearly, we show 
in fig. El the magnetic gluon propagator for the lowest 
two Matsubara modes ng = 0, 1 as function of the spatial 
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Fig. 16. Gluon propagator at jj, — 0.7, j — 0.04, together with 
the fj. = j = propagator for comparison. 



monientuni |q| for different fi. We see tliat there is very lit- 
tle change up to the deconfinement transition, after which 
there is a dramatic infrared suppression and clear ultravi- 
olet enhancement of the propagator. These changes makes 
it look more like an ordinary massive boson propagator, 
i.e., the gluon has acquired a magnetic mass which grows 
as fi increases beyond /x^. This effect does depend on the 
diquark source j, becoming weaker as j — > 0, as figure [T^ 
shows. Simulations at smaller j in the deconfined phase 
will be necessary to firm up the picture, but it seems un- 
likely that the j dependence seen at /i = 0.7 will be suf- 
ficient to wholly cancel the infrared screening effect. No 
dependence on j is seen at ^ = 0.3, 0.5. 



Figure ^1 shows the electric propagator for the low- 
est two Matsubara modes, for selected values of /i. While 
the electric propagator, as seen previously, remains un- 
changed up to /i = 0.7, it too is clearly suppressed in the 
infrared above deconfinement (although to a smaller ex- 
tent than the magnetic case), but appears unchanged in 
the ultraviolet. No dependence on the diquark source is 
seen. The same effects can be seen in both static (go = 0) 
and non-static propagators, although for no > 2 the split- 
ting between electric and magnetic sectors becomes harder 
to detect. The case of the static magnetic propagator is 
particularly interesting, since this is unscreened to all or- 
ders in perturbation theory, yet we observe a clear screen- 
ing effect. The increase of magnetic screening with j also 
suggests a non-perturbative origin, hinting at a relation to 
the non-vanishing (qq) condensate, rather than the pres- 
ence of light quasiquark degrees of freedom at the Fermi 
surface. The data in this case may be distorted by possi- 
ble finite volume effects, although the consistency between 
the hq = and hq — 1 data indicates that these do not 
affect the qualitative picture. 
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Fig. 17. Magnetic gluon propagator for qo = (top) and 
aqo = 2sin7r/16 (bottom) and diquark source j = 0.04, as 
a function of spatial momentum \q\ for various values of the 
chemical potential. 

5 Discussion and outlook 

Our results indicate that, at least for the (rather heavy) 
quark mass employed in this study, QC2D has three sep- 
arate phases: 

1. A vacuum phase, for fi < fio ~ "^tt/S, where the 
baryon density remains zero and all other physical 
quantities (with the likely exception of the hadron spec- 
trum [22) are unchanged. 

2. A confined, bosonic superfluid phase, for ^o < ^J■ < IJ-d-, 
characterised by Bose-Einstein condensation of scalar 
diquarks. The thermodynamics of this phase can be 
qualitatively (but not quantitatively) described by chi- 
ral perturbation theory, and the static quark potential 
is screened by the dense medium. 

3. A deconfined phase, for fi > fj,d, where quarks and 
gluons are the dominant degrees of freedom. In this 
phase, a Fermi surface of quarks is built up, leading to 
scaling of thermodynamic quantities of the same form 
as predicted by the Stefan-Boltzmann form, but with a 
different constant of proportionality. We interpret the 
latter as evidence for a non-zero binding energy in this 
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Fig. 19. Magnetic gluon propagator at jj, — 0.7 for various 
values of the diquark source j. 



phase, as kp > Ep. We observe Debye screening of the 
electric gluon propagator, as well as strong screening 
of both static and non-static magnetic gluon modes. 

The deconfinement transition occurs at apid ~ 0.65, which 
in physical units corresponds to /id « 600 MeV. The corre- 
sponding quark density may be estimated as Uq w llfm^ 



T" 



, e /« raw data 

, e /n corrected data 
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+'H 
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Fig. 20. The ratio Eq/uq vs. ^l, using both the raw data of 
table 01 and rescaled by a factor n^sB/^^sk 



directly from table^J or from fig. where lattice artefacts 
are taken into account, as «6.5fm^'^. 

Let us discuss some phenomenological implications of 
our results. First consider the nature of the deconfinement 
transition. Since both the low density "nuclear matter" 
and high density "quark matter" phases are characterised 
by a U(1)b breaking superfluid condensate, it is tempt- 
ing to postulate the quark/hadron continuity proposed for 
QCD with Nf = 3 light flavors in "30". In that case, the 
spectrum of physical excitations in both nuclear and quark 
phases is qualitatively similar and can be matched us- 
ing phenomenological insight, for instance that the phys- 
ical state with the smallest energy per baryon is the 6- 
quark bosonic state known as the H-dibaryon, which can 
Bose-condense to form a superfluid at nuclear densities. In 
QC2D, however, no matching is possible. The spectrum in 
the confined phase consists entirely of bosonic states, be 
they qq mesons, qq baryons and their conjugates, or glue- 
balls. The deconfined phase by contrast has quasiquark ex- 
citations with half-integer spin and mass gap Z\, the same 
energy scale as all non-Goldstone mesonic and baryonic 
excitations. The only states which are conceivably lighter 
are spin-1 quasigluons, with no counterpart in the low- 
energy spectrum predicted by xPT 23 . Since the spec- 
trum is discontinuous, it seems natural to suggest that 
deconfinement in this case is a genuine phase transition, 
rather than a crossover. 

The most intriguing of our results concern the mag- 
netic gluon propagator shown in fig. 1171 The long-distance 
screening observed in the static limit gp "*■ is not pre- 
dicted at any order of perturbation theory. However, this 
breakdown of perturbation theory should not be unex- 
pected as the magnetic gluon is an intrinsically non-per- 
turbative object. We have no explanation for the screening 
effect, but note that absence of magnetic screening in the 
static limit is a crucial ingredient of the celebrated calcu- 
lation of the gap scaling A oc exp(— 37r^/A/2g) predicted 
for high-density QCD [H]. 
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Finally, it is fascinating to speculate on the astrophys- 
ical consquences of our results. Figure EHl plots the energy 
per quark Eqjuq against /x, and shows a shallow minimum 
at \xa ~ 0.8. This feature occurs whether or not correc- 
tions for lattice artefacts are applied, and appears to be a 
robust prediction of our simulation. Inclusion of the gluon 
contribution Eg shown in fig. will ensure de/duq > as 
/z — > oo, as of course will recovery of the Stefan-Boltzmann 
scaling H22I23I) as asymptotic freedom sets in. Moreover, 
the minimum is not predicted by xPT, where e increases 
monotonically with Uq (see ea. (|2(J|) '). The minimum resem- 
bles a property known as saturation in nuclear physics, 
and implies that objects formed from a fixed number of 
baryons, such as stars, will assume their ground states 
when the majority of the material lies in its vicinity. Since 
the minimum lies above the deconfining transition, we de- 
duce that two-color stars are made of quark matter. 

Orthodox models of quark stars [32| are based on a 
simple equation of state such as the Bag Model, which 
predicts a sharp first order deconfining transition. The 
resulting stars have a sharp surface where p = 0, along 
which quark matter coexists with the vacuum. In QC2D 
by contrast, the state of minimum e/uq has p > 0. A 
two-color star must therefore have a thin surface layer, 
perhaps better described as an atmosphere, formed from 
scalar diquark baryons, and whose density falls continu- 
ously to zero as the surface is approached. Matter in the 
range 0.4 ^ ^a < 0.7, is mechanically unstable, since ac- 
cording to iig.\^dp/de < and the resulting sound speed 
imaginary. At the base of the diquark atmosphere there 
will be a sharp increase in both pressure and density, and 
the bulk of the star will be formed from quark matter 
with /x > /id. Precise radial profiles, and the relation be- 
tween the star mass M and its radius i?, must await more 
quantitative information of the equation of state, which 
requires the correct normalisation of Sq and Eg . 

We finish by outlining future directions of study. The 
lattice spacing used here is quite coarse, and as illustrated 
in section \T?B lattice artefacts are quite substantial. It will 
be important to repeat the study on finer lattices in or- 
der to gain control over these artefacts. This is currently 
underway. 

It will also be desirable to find the correct rescaling 
factors for the energy densities, via a non-perturbative 
determination of the Karsch coefficients using simulations 
on anisotropic lattices |33II34| . 

Since the differences between QCD and QC2D are great- 
est in the chiral limit, the heavy quark mass employed 
here can be considered a positive rather than a negative. 
Nonetheless, it would be desirable to study a system with 
lighter quarks, to explore the mass dependence of our re- 
sults and interpolate between the regime we have been 
exploring and the chiral regime. Better control over the 
limit j —f is also required. 

Beyond the issues considered in this paper, we are in- 
tending to study the hadron (meson and diquark) spec- 
trum and the fate of Goldstone as well as non-Goldstone 
modes in the dense medium. Issues of interest there are 
the fate of the vector meson and the possibility that it be- 



comes light in the dense phase |35ll3(illll| . Also of interest 
is the pseudoscalar diquark, which may provide a pointer 
to the restoration of the U(1)a symmetry in the medium 
pTfj . We also intend to study the Gor'kov quark propaga- 
tor in momentum space, which will provide information 
about effective quark masses and the superfluid diquark 
gap. Gauge-invariant approaches to identifying the pres- 
ence of a Fermi surface in a Euclidean simulation can also 
be employed ^Q. Finally, once ensembles of the dense 
medium are available on a fine lattice, it will be interest- 
ing to analyse topological excitations ^21i to see to what 
extent deconfinement in dense baryonic matter resembles 
deconfinement in a hot medium. 
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